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I. INTRODUCTION 



The asymmetric simple exclusion process (ASEP) is one of the fundamental models in non-equilibrium statistical 
mechanics. The ASEP involves particles that perform asymmetric jumps on a discrete lattice under the exclusion 
constraint: two particles cannot occupy the same site at the same time. This very simple and minimal system 
appears as a building block in more realistic descriptions for low-dimensional transport with constraints. Such phe- 
nomena occur in various contexts and scales ranging from micrometric cellular motors to traffic networks |32l 140) . 
The remarkable properties of the ASEP and its numerous variants have stimulated hundreds of studies during the 
last two decades. The ramifications of this mundane-looking model through non-equilibrium statistical mechanics, 
combinatorics, probability, random matrices and representation theory are tremendous |T1 151 1^ 15^ [5^ H5] . 

Generically, a system out of equilibrium carries at least one non- vanishing current in its steady state (this is related 
to the breaking of detailed balance). Such currents can be considered as archetypal observables for non-equilibrium 
behaviour j40j . Classifying the different independent stationary currents in a given system, identifying some generic 
features and relations obeyed by the distributions of these currents and calculating their statistical properties as 
functions of the control parameters, are some of the important tasks in non-equilibrium statistical physics. Current 
fluctuations are usually non-Gaussian: their characteristics can be quantifled by the moments of the current (mean- 
value, variance, skewness, kurtosis...) or by a large-deviation function that measures the probability for the current 
to assume a non-typical value. There is growing evidence that large-deviation functions play a crucial role in non- 
equilibrium statistical physics, akin to that of thermodynamic potentials at equilibrium [101 144j . 

Currents transport information from one part of the system to another. In particular, in a system far from 
equilibrium, boundary conditions can drastically alter the behaviour of the bulk even if interactions are short-ranged 
(in contrast with the generic case at equilibrium). For example, as was recognized in the earliest studies [551 HH]: phase- 
transitions can be induced by the boundaries or by a localized alteration of the dynamical rules. It is therefore crucial 
to specify the boundary conditions (e.g. periodic, twisted, open boundaries, inflnite system...). The phenomenology 
of the system and the mathematical techniques that are used to analyse it depend strongly on the chosen boundary 
conditions. 

In the present work we consider the exclusion process on a finite lattice with open boundaries, which can be viewed 
as a model for a conducting rod in contact with two reservoirs that are not in thermodynamic equilibrium with each 
other (for example, they are at different temperatures, or have different chemical, mechanical or electrical potentials). 
Besides, an external field may be applied to the system. This external field and the reservoirs drive a current through 
the system and, in the long time limit, the connecting rod reaches a non-equilibrium stationary state. Our aim is 
to study the statistics of this stationary current. We give explicit formulae, equations (151 to (19), for the cumulant 
generating function of the totally asymmetric exclusion process (TASEP) with open boundaries that are valid for 
arbitrary values of the entrance and exit rates and for all values of the system size L. We emphasize the fact that 
our results are of combinatorial nature and not only asymptotic: they describe the TASEP in all possible regimes, 
including the phase-transition lines. This is in contrast with the asymptotic expression for the large-deviation function 
of the current obtained very recently using the Bethe Ansatz [H]. Indeed, as of today, the Bethe Ansatz for the open 
TASEP [7] seems to be tractable only in the L — )• oo limit, and only inside the low and high density phases far from 
the phase-transition lines. We show that the formula obtained in [8] for L -> cx) can be retrieved as a limiting case of 
our general results. The parametric representation we have found is obtained by using the Matrix Ansatz technique 
[TJ [TS] . The calculations are very cumbersome and some combinatorial patterns have been guessed rather than fully 



2 



calculated. The formulae we obtain must therefore be considered as conjectures. But as we shall explain, the results 
are exact albeit they have been obtained, at present, in a non-rigorous manner. They have also been thoroughly 
verified against exact computations for systems of small sizes. Besides, all known special cases can be deduced from 
our general result. 

The outline of this work is as follows. In section [llj we recall the definition of the model and the basic properties 
that will be used later (the master equation, the matrix solution and the phase diagram). Section III contains the 
main results of this work; we first restate the problem of current fluctuations in the standard mathematical framework, 
and then describe the simpler case where all the parameters in the system are equal and taken to be 1. The general 
formulae are presented in IIIC and are followed by a physical discussion of the behaviour of the system through 
the phase diagram; we also discuss some connections with previously known results. In section |IV[ we explain the 
line of reasoning that we have followed to obtain the formula for the current fluctuation, show where the unproven 
assumptions reside, and describe some numerical veriflcations. The last section is devoted to concluding remarks. A 
detailed derivation of the large L limit [8 from the general formula is given in the appendix. 



II. DEFINITION OF THE TASEP AND BASIC PROPERTIES 



We consider the totally asymmetric simple exclusion process (TASEP) on a flnite lattice of size L with open 
boundaries. Each site of the system can be occupied by at most one particle [exclusion condition). The dynamics of 
the model is defined by the following stochastic rules (see figure [T]): a particle at a site i in the bulk of the system 
(with 1 < i < i — 1) can jump with rate p — 1 (i.e with probability dt during the time interval dt) to the site i + 1 if 
this target site is vacant; if site 1 is empty, a particle can enter with rate a; a particle at site L can leave the system 
with rate f3. The entrance and exit rates represent the coupling of the finite system with infinite reservoirs located 
at its boundaries. At a given time, the system is in one of its 2^ possible configurations and evolves according to its 
stochastic dynamics as a Markov process. The evolution of the system can be encoded in the Markov matrix M as 
follows: the probability Pt{C) of being in configuration C at time t satisfies the master equation 

^=^M(C,C')F.(C'). (1) 

The non-diagonal matrix element Af(C,C') represents the transition rate from C to C. The diagonal part M{C,C) = 
— X^c'^C M{C' ,C) represents (minus) the exit rate from C. The Markov matrix is a stochastic matrix: the sum of the 
elements in any given column vanishes. 




FIG. 1: Illustration of the TASEP with open boundaries on a finite lattice with L sites. 



In the long time limit, the system reaches a steady state in which each of the 2^ possible configurations occurs 
with a stationary probability. This steady state probability lies in the kernel of the Markov Matrix: the rules of 
the ASEP ensure that this kernel is non-degenerate and that all other eigenvalues of M have strictly negative real- 
parts that correspond to relaxation states with a possible oscillatory behaviour (Perron- Frobenius theorem). Finding 
this stationary measure is a non-trivial task: the model is far from equilibrium with a non-vanishing steady-state 
current, there is no underlying Hamiltonian and no temperature. Therefore the fundamental principles of equilibrium 
statistical mechanics, such as the Boltzmann-Gibbs law, cannot be used. 

The exact calculation of the stationary measure for the TASEP with open boundaries and the derivation of its phase 
diagram have played a seminal role by triggering a whole field of research on exactly solvable models in non-equilibrium 
statistical mechanics. We recall that the fundamental observation [T^ is the existence of recursion relations for the 
stationary probabilities between systems of different sizes. These recursions are particularly striking when a — (3 = 1 
p3j : they can be generalized to arbitrary values of a and /3 [ini 112] and also to the more general case in which 
backward jumps are allowed (PASEP). The most elegant and efficient way to encode these recursions is to use the 
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Matrix Ansatz [TB]. A configuration C can be represented by the binary string of length L, (ti, . . . , tl), where = 1 
if the site i is occupied and = otherwise. Then, to each C, the following matrix element is associated: 

m-^(«in(T.D + (l-r,)E)|^). (2) 



i=l 



The scalar P{C), thus defined, will be equal to the stationary probability of C if the operators D and E, the bra- vector 
(a I and the ket-vector satisfy the following algebraic relations 



DE = D + E 

1 , 



HE = i(a|. (3) 

a 

This algebra allows us to calculate any matrix element of type ([2|. The normalisation constant in equation ([2| is 
given by 

ZL^{a\iB + Ef\(5). (4) 

For a = /? = 1, is a Catalan number [IB]. More generally, the Matrix Product Representation method has proved 
to be very fruitful for solving many one-dimensional systems: a very thorough review of this method can be found in 

m- 

From the exact solution, the phase diagram of the TASEP, as well as stationary equal-time correlations and density 
profiles, can be determined. In the limit of large system sizes, the phase diagram (figure [2]) consists of three main 
regions 

• For a < min(/3, 1/2), the system is in the Low-Density Phase and its behaviour is driven by the entrance rate 
a. The bulk-density is p = a and the average current J = a(l — a). 

• The High Density phase, for /3 < min(Q!, 1/2) is characterized by p = 1 — /3 and J = /3(1 — /3). 

• In the Maximal Current Phase, with a > 1/2 and /3 > 1/2, the bulk behaviour is independent of the boundary 
conditions and one has p = 1/2 and J = 1/4. However, correlations with the boundaries decay only algebraically 
whereas they decay exponentially in the two other phases. 

• The low and high density phases are separated by the 'shock-line', a = /? < 1/2, across which the bulk-density 
is discontinuous. In fact, the profile on this line is a mixed-state of shock-profiles interpolating between the 
lower density pi = a and the higher density pn = 1-/3. 

Detailed properties of the phase diagram are reviewed for example in [U [HI HI] . We note that the phase diagram 
was obtained in [28] through physical reasoning by using a hydrodynamic limit and mean-field arguments (see also 
[43]). However, a finer analysis does require the knowledge of the exact solution. 



III. CURRENT FLUCTUATIONS FOR THE TASEP 



Many properties of the ASEP have been understood using different techniques (Matrix Ansatz, Bethe Ansatz, 
Random Matrices) [J [HI [211 [HI 1321 [55WT1 1^ . However, the determination of current fluctuations in the original 
TASEP model with open boundaries has remained a vexing and challenging unsolved problem. This question is 
interesting and important: first, in presence of reservoirs, the model is quite realistic and can be related to real 
experimental situations [321 [57] (a niore detailed discussion can be found in second, the exact calculation of 

fluctuations to all order is akin to determining the large-deviations of the current which are expected to play a central 
role in non-equilibrium statistical mechanics jlOl I44j. 



A. Statement of the problem 



We consider the TASEP with open boundaries and want to study the total (i.e., time-integrated) current that has 
flown through it in the long time limit. One way to quantify this total current is to place a counting variable Nt at 
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FIG. 2: The phase diagram of the TASEP as a function of the boundary rates. 



the entrance site. At t — we have A'o — 0; each time a particle enters the system, we increment the value of Nt by 
1. Hence, Nt is a random variable that counts the total number of particles that have entered the TASEP between 
time and t. Because the size L is bounded and no particles are created or destroyed in the bulk, Nt also represents, 
when < — >■ oo, the number of particles that have crossed any bond in the system, or have exited from the TASEP from 
its right boundary. We call Nt the total current at time t and we intend to study its statistical properties. 
When t — >■ oo, the expectation value of Nt/t converges towards the average stationary current J: 

lim^ = J(«,/3,L) = %^. (5) 

The value of Zl, defined in equation Q, was determined exactly in [TB] using the Matrix Ansatz. 

The variance of Nt in the long time limit also increases linearly with time. This allows to define a 'diffusion-constant' 
A, as follows: 

limi^^^l^M!.A(«,AL). (6) 

t->-oo t 

An exact expression for A(a,/3, L) was obtained in |17j . It involved an extension of the matrix method to represent 
some two-time correlations in the stationary state. For a — f3 — I, the formula for A is quite elegant and involves 
simple factorial factors (see the next section) . Unfortunately, the general result is somehow less compact (see equations 
(58) to (61) in [H]). 

More generally, one can define moments and cumulants to all orders for the random variable Nt- This data can be 
encoded in the exponential generating function (cxp(7 A^t)). The moments of Nt are obtained by taking successively 
the derivatives of this function at 7 = 0. In the long-time limit, we have 

{eMlNt))c^eMEh)t). (7) 

(A mathematical equality is obtained by taking the logarithm of both sides, dividing by t and taking t — >■ 00). The 



derivation of this result is recalled in Section IV Because the logarithm of (exp(7 A^t)) generates the cumulants of Nt, 
we observe that these cumulants grow linearly with time and that their values are given by the derivatives of E{'j) 
at 7 = 0. The cumulant generating function £'(7) is related to the large-deviation function of the current by Laplace 
transform. For finite-size systems, both functions carry the same information. 

The function E{j) was calculated for the symmetric exclusion process in [M]. Very recently, the analysis of the 
Bethe Ansatz equations, for L — >■ 00 was carried out for the asymmetric case, in the low and in the large density 
phases [5] . In the present work, we obtain an explicit representation of the generating function of the cumulants of 
the current for all values of a and /3 and for all values of the system size L. The technique used is an extension of the 



matrix method (see section IV ) 



For pedagogical reasons we shall first discuss the case a — (3 — 1, which belongs to the maximal current phase. 
Here, the formulae are explicit, quite simple and appealing. The general case with arbitrary values of (a,/3) will be 
presented in a separate subsection. 
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B. The cumulant generating function for a — /3 = 1 



Here, we consider the values a — /3 — I. Historically, the TASEP with a ~ /S ~ I was the first case for which the 
stationary measure was determined exactly |13] . For the cumulants of the current, these special values lead to simple 
mathematical expressions. 

The generating function of the cumulants E{-^) is given as a function of 7 through the following representation in 
terms of a parameter B: 



E 

fc=i 



(2fc)! [2k{L + l)]\ 
kl [fc(L + l)]![fc(L + 2)]! 2fc ' 

(2fc)! [2fc(L + l)-2]! B'' 



(8) 



kl [fc(i + l)-l]![fc(L + 2)-l]! 2fc ■ ^' 

By expressing B in terms of 7 in equation ([s]) and substituting in equation (|9| we can calculate as a function of 7 
to any desired order. The coefficients of this expansion that we denote by 

ii;(7) = i?i7 + 7' + Itt' + ■ • ■ , , (10) 

give the successive cumulants of Nt in the long time limit. For example, we obtain 

limM^j^^,^^ + 2 

t-s-oo t 2(2L + 1) ^ ' 

This expression is identical to that found in [16 . When L — > 00, J — >■ 1/4 as expected in the maximal current phase. 
At the second order, we obtain 

(N?) - (N,)^ 3 (4£ + l)![L!(L + 2)!]2 

}^ -t =^ = ^2=2 [(2L + l)!]3(2L + 3)! ' ^^^^ 



This is the same formula as the one derived in Jlj- When L 00, we have A ~ ^ gl^ L"-'-/^ i.e., the diffusion constant 
vanishes in the maximal current phase. 

The third cumulant, known as the skewness, is given by: 

^ [(£ + l)!]2[(£ + 2)!]4 r (L + l)!(L + 2)!(4L + 2)!(4£ + 4)! (6L + 4)! ^ 

^ (2L + l)[(2L + 2)!]3 I (2L+l)![(2L + 2)!]2[(2L + 4)!]2 (3L + 2)!(3L + 6)! J ' ^ ' 

This formula was not known before. For a large system, the skewness behaves as 

S3.HH__J|2^__0.0090<,78... (14) 

Carrying on the elimination of B, the next few orders can also be found explicitly. It is found that the k-th cumulant 
scales as 7r(7rL)'^/2-3/2 Jq^. > 2. We determined the constant prefactor for the first few values of k. 

The fact that the large deviation function is obtained in parametric form should not come as a surprise. On the 
contrary, this structure seems to be rather common; to our knowledge it appeared first in |18j . where the exact 
large deviation function for the TASEP on a ring was calculated by Bethe Ansatz; it also occurred in the case of a 
defect particle on a ring [TS]. The complete solution for the current fluctuations of the ASEP on a ring involves a 
tree-structure that is again written parametrically [35 • More generally, it was shown in [31 131 [TUI that for a large class 
of non-equilibrium diffusive systems, characterized by a linear conductivity and equilibrium diffusion coefficient, the 
large deviation function can be expressed through a parametric set of equations. 



We observe that the limiting value of the skewness, given in equation (14), is a finite number. More generally, one 
can construct simple ratios of cumulants that have a finite limit when L — > 00. It should be possible to 'measure' 
such numbers in simulations and to test the universality of these ratios in a manner similar to |12j . 
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C. The general case 



We now give the result for the cumulant generating function, vahd for any system size L and for arbitrary values 
of a and 13. The mathematical structure is the same as above: 



7 



oo 



2k ' 



(15) 



(16) 



fe=i 



The coefficients Ck and of the series are functions of the rates a and /3 and of L. Their explicit expressions are 
given by 



Cfe(a,/3) 



{O.a.b} 



dz F{zY 
2i7r z 



r dz F(z)'' 
and Dk{a,p) = f — 

J{o,a.b} 2«7r (1 + z)2 



where the rational function F(z) is given by 

-(1 + Z)2i(l-Z2)2 



Fiz) 



z^{l - az){z - a)(l - hz){z - b) 



with a 



1 - a 



and h 



1-/3 

/3 



(17) 
(18) 

(19) 



Equations ( 15 ) to ( 19 ) provide an exact representation for the generating function of the cumulants of the current 



in the TASEP with open boundaries. These equations are valid for any system size L and for arbitrary values of 
the parameters a and (3. Before we proceed further, we explain the meaning of the symbol ^^^y. it represents the 
complex integral along 3 infinitesimal contours that encircle the points 0, a and b in the complex plane. Equivalently, 
using the Cauchy formula we could have written equations (171 and (18) as 



C7fc(a,/?) 
Dk{a,p) 



Residue 
Residue 

20— 0,a,b 



F{zf 



Zq 



F{z)^ 



We note that similar integrals have already appeared in closely related problems [6lll5j. Here, the notations implicitly 
assume that 0,a and h are distinct. If two or three of these numbers coincide, one should take into account the 
corresponding residue only once. Hence, for a = /3 = 1, which implies a = 6 = 0, only the residue at must be 
calculated. The resulting explicit expressions for Cfc(l, 1) and Dk(l, 1) give the parametric representations ^ and ^ 
that were discussed in the previous section. 



By inverting the series (15 1 order by order and substituting the result in (161, we can derive closed expressions for 



the first few cumulants of the current. For example, the mean- value of the current is given by 



J 



Ci{a,P) 



(20) 



This formula is, of course, the same as obtained in We emphasize that Ci(a,/3) and Di{a,f3) coincide with 

Zl and respectively, as can be seen by comparing equations (171 and (18) for fc = 1 with equation (BIO) of 

Reference [16] (see remark [45]). 

At second order we obtain an expression for the diffusion constant 



A 



CI 



(21) 



We remark that C2 and D2 are natural generalisations of Ci and Di. In this form, the diffusion constant looks more 
compact than the formula found in [T7] . The two expressions must coincide but we have not endeavored to prove this 
fact analytically for all values of a,/3 and L. We can carry on this procedure further to a few more orders to obtain 
higher cumulants and we did so with the help of a symbolic mathematical calculation tool. 

It is of greater interest to analyse the behaviour of the large-deviation function in the different phases of the model 
when L — ^ 00. 
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In the Maximal Current phase, with a > 1/2 and (3 > 1/2, the parameters a and b lie inside the unit circle. 
Hence, the contour integrals that appear in equations ( 17 1 and ( 18 ) can be replaced by a single integral along the 



unit circle. Then, we apply the saddle-point method to estimate the asymptotic behaviour of Ck and when 
L — >■ oo. The saddle point is at z = 1 and we observe that the values of a and b do not influence the saddle-point 
estimation at dominant order: in fact they contribute by the same multiplicative factor [(1 — a) (1 — 6)]^*^, that 
can be reabsorbed in the parameter B. Therefore, the behaviour of the cumulants in the large L limit does not 
depend on the boundary rates a and /3 in the Maximal Current phase (as expected). The results at dominant 
order are the same as those obtained for a = /3 = 1, the special case discussed in the previous section. 

• In the Low Density phase, a < min(/3, 1/2), the parameter a is outside the unit circle, we have a > b and the 
position of b with respect to the unit circle is not determined. In the large L limit, it is the pole at z = a that 
contributes dominantly to the values of Ck and Dk and the parametric representation (15) and (16) becomes 
(see the appendix for more details) 



1 d 



E 



k=l 

^ i-! dzfe-U (z + 1)2 i 



k=l 



where the function (/'(z) is given by 

<j>{z) = (z - a)F(z) 



2\2 



-(l-fz)2^(l-z2) 
z^{l - az){l - bz){z- 



(22) 
(23) 

(24) 



These expressions can be used to calculate the first few cumulants in the low density phase: 



El 


= P(l 


-P) 


E2 


= P(l 




E3 


= P(l 






= P(l 




E5 


- p(l 




= a. 







240p^ + I20p'^) etc. 



(25) 



with the mean density p 

In fact, using equations (22) and (23), the function for E{'^) can be obtained in a closed form thanks to the 
Lagrange Inversion Formula |22| , as explained in the appendix. This leads to 



Eh) 



1 et 



(26) 



This expression is identical in the TASEP case to the one obtained by Bethe Ansatz in 8]. (In [8 , the general 
PASEP is studied: this adds a prefactor [p — q) to (26), where p and q are the rates of forward and backward 
jumps, respectively, and modifies the definition of a). We remark that the limit formula (26) is rather simple 
and that it is totally independent of the specific form of the function (j){z), as can be seen from the derivation 
given in the appendix. An elementary and physical derivation of this result can be given by using macroscopic 
fluctuation theory [U [11] . 

• The High Density phase is symmetrical to the low density phase under the exchange of a and /?. Therefore a 
separate discussion is not required. 

• The shock line a — {3 < 1/2, i.e., a^b>l can also be analyzed from our general formula. We have calculated 
explicitly the first few cumulants and obtained the following scaling behaviour 



Ek tka{l - a){l - 2af 



^ for fc > 2 . 



We recall that the current is given by Ei = a{l — a). The numerical coefficients are given by £2 = 2/3, 
£3 = —1/30, £4 — 2/315, £5 = —1/1890... The fact that the higher cumulants grow without bounds with L 
whereas they are bounded in low and high density phases may come as a surprise. However, it is known that 
in the shock phase, particles have a vanishing chemical potential and the equivalence between canonical and 
grand-canonical descriptions breaks down [551 ISS] • 
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More precisely, this puzzling scaling can be understood from the domain-wall picture, introduced in [T7] to 
describe the discontinuity by a factor 2/3 in the diffusion-constant along the shock line. In the limiting case 
a — /3 1, the dominant configurations are shock profiles between a low density region with pi — a and a high 
density region with ph = 1 — (3. The shock is localized on a single site. The dynamics becomes equivalent to that 
of a symmetric random walker on a lattice oi L + 1 sites and confined by two reflecting boundaries. The number 
of particles having entered the TASEP corresponds to the number of leftward steps of the shock. The statistics 
of the steps performed by an effective random walker between two reflecting walls is a well-posed problem that 
can be studied independently. For instance, we find that the skewness (third cumulant) of the random walker 
grows as — L/30. However, if we introduce some bias in the jumping rates (which corresponds to the high or low 
density phase) the walker becomes localized near one of the walls and all cumulants remain finite when L — > (X). 
Finally, the e^'s can be calculated by remarking that the random walker between two reflecting walls can be 
mapped to the TASEP with 2 identical particles on a periodic lattice, and then using the results of [T5]. 

IV. OUTLINE OF THE ANALYTICAL PROCEDURE 



We now describe briefly how the expressions ( 15 ) to ( 19 1 were obtained. Detailed explanations are deferred to a 
forthcoming article. 

A. General setup 

In order to study the statistics of Nt, the total number of particles that have entered into the system, one introduces 
the joint distribution Pt(C, N), the probability of being at time t in configuration C and having Nt = N . The evolution 
equation of Pt{C, N) can be written using the Markov equation Q. It is then useful to introduce the Laplace transform 

F*(C)=^c''^Pt(C,7V). 
iv 

These generalized weights satisfy a deformed Markov equation 

^=M{^)F,. (27) 

The matrix Af(7) of size 2^ is given by 

Af(7) = M + (e'^ - l)Mi , (28) 

where M is the original Markov operator and the matrix Mi contains only those transitions in which a particle 
enters the system, i.e., Afi (C,C') = a if C evolves into C by adding a particle at site 1 and Mi(C,C') ~ otherwise. 



Equation (27) is formally solved as Ft = exp(M(7)t)Po which, in the long time limit, behaves as 

Ft e^(^)* \E{^)) , 

where -£(7) is the dominant eigenvalue of ^(7) and |-B(7)) the associated eigenvector (it is unique for sufficiently 
small 7). Thus, we obtain equation ([t]) 

(exp(7iVt)) = - eME{l)t) . 

c 

The cumulant generating function is therefore identical to the largest eigenvalue of the deformed operator M{'~f) and 
the problem of determining the statistics of Nt has been traded for a spectral problem. This question can be tackled 
using different methods. Bethe Ansatz is one technique for integrable systems. Another approach, valid for small 
values of 7, is to perform a perturbative expansion around the dominant eigenvector, |0), and the dominant eigenvalue, 
0, of the original Markov matrix M: 

E{^) = S17 + 1^7' + 1^7' + ■ • ■ 
\E{l)) = |0)+7|1)+7'|2) + ... (29) 
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Carrying out the perturbative expansion explicitly, we find that the fc-th order correction to the eigenvector satisfies 
an equation of the type 

M|A:) =7^fc(|0),...|A:-l)) , (30) 

TZk being a linear functional. For example, we have 

Af|0) = , M|l) = {El - Afi)|0) etc ... 

Moreover, the solvability condition (obtained by using the fact that vector (0| = (1,1,...,!) is the left null-eigenvector 
of the Markov matrix M) allows to express the k-th term in £^(7) (i.e., the cumulant of order k) as a linear function 
of the k — 1 vectors |1), |2), . . . |fc — 1). For example, we have 

Ei^{0\Mi\0), ^2-2(0|Mi|l) + (l-2Si)(0|l) etc... (31) 



The cumulants of the current can be found if we are able to solve the set of linear equations ( 30 1 for all values of 
A: > 0. At order 0, the stationary state |0) and the mean current Ei were calculated thanks to the Matrix Ansatz, 
recalled in equations ^ (|3| and Q. In [T7], the first order correction was also obtained by a Matrix Ansatz: the 
required algebra was constructed by taking tensor products of two quadratic algebras of the type ([3]). This allowed 
to calculate the diffusion constant E2- 

B. Steps of the calculations and Numerical tests 

The computations that we have carried out can be summarized by the following steps: 

(i) The Matrix Ansatz used in |17] has been simplified and generalized to all orders. The operators required to 
calculate the (fc + l)-th cumulant are denoted by and E^. These operators are constructed by using 2fc + 1 tensor 
products of the original D and E's. This may appear daunting at first sight but in fact this is not so bad: the Dj, 
and Efc already appear in the studies of multispecies exclusion processes [101 1301 El] ■ There they were introduced as 
formal objects to construct the stationary measure. Here, they are used as tools for calculations. We have also found 
the boundary vectors (ak| and |/3k)- 



(ii) The fact that and E/j allow to solve the system (301 at each order has been proved recursively. Using the 



relations (31 1, an expression for the cumulant Ek in terms of matrix elements involving and E^ and the previously 



determined Ej^s for j < k can be written. 

(iii) The next step is to calculate the matrix elements. Typically, one has to determine a 'normalization'-type term 
(afe|(Dfc + Ek)^\Pk) that generalizes Such a matrix element can be found by using the image method [T71[Tn] in 
a space of dimension 2k + 1. The resulting expression involves sums of products of binomial factors. These binomial 
factors can be expressed as multiple contour integrals and the total sum can be recast as a determinant. As this stage, 
the cumulants are written as complex integrals over determinantal expressions. 

(iv) The diffusion constant {k = 2) was recalculated. The skewness {k = 3) was then determined by evaluating the 
contour integrals. These integrals produce a 'generic' term and a very large number of boundary terms. It was found 
that the global contribution of the boundary terms cancels out for A; = 2 and 3. It was also observed that a parametric 
expression for 7 and E at order 3 as a function of an arbitrary parameter B allows to retrieve the diffusion constant 
and the skewness. 

(v) At k-th order, a calculation of the 'generic' term was performed by assuming that the various contributions of 
boundary terms globally cancel out. It was observed that the resulting formula could also be obtained through the 



parametric expressions (15) and (16 1 for 7 and E at order k with respect to B, that were guessed at all orders. 



We admit that our final result was guessed rather than fully worked out. The gaps in our derivations (steps (iv) 
and (v)) have to be filled. However, we emphasize that we have been using a systematic procedure and that we are 



absolutely sure that the equations (15) to (19) are correct. We have tested numerically the conjectured formulae in 
the following cases: 
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For a = /3 = 1, the formula (131 for the skewness was compared with the exact result obtained by solving the 
linear equations (30) for the first three orders for systems of sizes L < 10. The results agreed in all cases. The 



expressions are rational numbers which involve integers having more than 20 digits. 
For a — /3 — 1, the first 6 cumulants were verified for L < 10. 



The formula for the second cumulant (21 ) as a function of a and /3 was tested against the exact result for L = 4. 



We have chosen arbitrary values of a and /3 throughout the phase diagram (a dozen of different cases). By 
taking rational values for a and /3, we insure that the formal mathematical program (Mathematica) performs 
exact calculations on integers when solving the linear equations ( 30 1 . We have worked with systems of size 



L < 10 and tested the formula up to the 6th cumulant. The integers that appear involve hundreds of digits (450 
in the worst chosen case): the expressions (151 to (19) give the correct answer in all cases. 



Finally, we emphasize the formulae given here allow us to re-derive results that were known previously, such as 
equations (12) and (26). This provides an indirect verification of their correctness. 



V. DISCUSSION 



In the present work, we give closed formulae for the generating function of the cumulants of the current in the 
open TASEP. These results hold for all values of the boundary parameters a and f3 and all values of the system size 
L. We emphasize the fact that our expressions are exact and not only asymptotic: they are of combinatorial nature. 
They allow us to describe the ASEP at all points of its phase diagram, including the phase-transition lines. The 
cumulant generating function is given in the form of a parametric representation - equations (15) to (19): a similar 
mathematical structure can be found in other works (THl UHl 131] and it can be related physically to the additivity 
principle [51 [TU]. 

The statistics of the current in the open TASEP had remained a challenging open problem for many years and no 
exact solution for the full distribution of the current was known for finite size systems. The Bethe Ansatz equations 
for the open ASEP were studied in but they are valid only on some surfaces in the parameter space: this restriction 
seemed to be a major obstruction to the computation of the large-deviation function. Only recently, a very subtle 
analysis of the Bethe equations, valid in the L — )■ cx) limit, together with some conjectures on the asymptotic locations 
of the Bethe roots, was carried out in [8], leading to an expression for the cumulant generating function. However, 
the result in [5] can only be established deep inside the low and the high density phases and it is an asymptotic 
expression, valid only when L — >■ oo. In our work, we have followed a different path and used the Matrix Product 
Representation ^ 116) to calculate the cumulants of the current order by order. We have performed some explicit 
calculations and uncovered the general structure of the solution. From the mathematical point of view, our formulae 
are only a conjecture but we have verified it in dozens of cases and derived from it all the previously known results. 
We have absolutely no doubt that our expressions are true. In particular, the main expression obtained in 8 can be 
derived as a limiting case of our results. 

We have decided to present the final formulae before completing the proof, because we find the results elegant and 
sufficient by themselves. Furthermore, they allow to draw some interesting physical consequences and to open new 
problems. We think that the proof is a question of carrying out a very long computation rather than having some 
deep mathematical insight. We are presently working on this aspect. We hope to have given enough details to the 
reader to clarify what kind of assumptions were made, that allowed us to jump to the final result and to guess the full 
structure of the solution. Another possibility, now that the final formula is known, is to search for a direct method to 
check it: after all, the cumulant generating function is nothing but the largest eigenvalue of a known operator. 



Besides completing the derivation, we also intend to extract from equations ( 15 ) to ( 19 ) the scaling limit of the large 



deviation function in the L ^ oo limit. It would be interesting to compare that scaling form with recent numerical 
results obtained by Monte-Carlo and DMRG methods 25, 31, 33j and also to investigate the crossover with theoretical 
results derived on the infinite lattice [5TJ [31] . 

Finally, we have considered here only the TASEP. But the matrix method can also be applied to the partially 
asymmetric case (PASEP), and we have checked that tensor products of the PASEP quadratic algebra that appeared 
in multispecies PASEP models |35j allow us to solve the hierarchy of equations (30) for the cumulants. We believe 
that the parametric representation still holds in the PASEP case and that combinatorial tree-structures akin to those 
found for the PASEP on a periodic ring in [34] will probably play an important role. Besides, complex integral 
representations for matrix elements analogous to those we have used here also appear in the PASEP [TJ ^ [371 ISH] ■ 
There must exist a general and hopefully elegant structure that encompasses all the cases, though we are well aware 
that such a structure may be difficult to discover. 
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Appendix A: Derivation of equation (26) 



Let Co,Ca and Cb be three infinitesimal contours (circles) that surround the points 0, a and h in the complex plane. 
According to equations (171 and (18 1, we must calculate contour integrals along these three circles. We denote by Fi 
the unit circle. 

In the low density phase a < min(;5, 1/2), the parameter a is outside the unit disk and we have a > b. The point 
b can be either inside or outside the unit disk. In order to perform a saddle-point analysis in the large L limit, a 
contour must pass through the point z — 1. We therefore deform the infinitesimal contour Co encircling into the 
unit circle Ti. But then we must remove the contribution of the pole 1/a of F{z) which is inside Fi (see figure |3|. 
We must also distinguish the cases 6 < 1 or 6 > 1: For & < 1 we can write, formally: Co + Ca + Cb — Ti — C-^ja + C'a- 
For 6 > 1 we can write Co + Cq + Cb = Fi — Ci/o — Ci/b + Ca + Cb- Now the function F satisfies: F{z) — F(l/z). 
Thus for zq = ^5 we have 



Residue 
Residue 



F{zf 



F{z) 



(z + iy 



— —Residue 
= —Residue 



F(z)'= 



, 1/2^0 



F{z)'' 



1/^0 



(Al) 



Formally this can be written as Ci/a 
to Fi + 2Ca if 6 < 1 and to Fi + 2Ca 



-Cq and Ci/f, 
2Cfc if 6 > 1. 



'Cb- Finally the contours for the complex integrals reduce 



1/a 



1/b 




1/b 







FIG. 3; The two sets of equivalent contours in the complex plane. In this figure, we have supposed that b lies inside the unit 
disk. 



The integral over Fi is estimated by saddle point: it is of order 4^. The integral over Ca is of order {2+a+l / a)^ 
(because a > 1). The residue at b contributes only when 6 > 1 but in that case it is of order (2 + 6 + l/b)^ <C 
(2 + a + 1/a)^ (because b < a). Hence, in all cases, in the large L limit, the dominant contribution comes from the 
contour around a, which has to be counted twice. Finally, using Cauchy's formula, we can write 



where the function ^(-z) has been defined in \2A\- Similarly, we have 



(A2) 



(A3) 
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Substituting these expressions into the formulae (15) and (16) leads to equations (22) and (23). 

We now need the Lagrange Inversion formula |22j, that can be stated as follows. Suppose that the two complex 
variables z and w are related in terms of the parameter B as 

(A4) 



(A5) 



w — z + B(p(w) 

where (j) is locally analytic. Then, any function G{w) can be expanded as a power-series in B as follows: 
GH = G(z) + ^0(z)G'(z) + — - {(0(z)f G'(z)} + . . . + (^^^ {{mf^^ G'(z)} + . . . 



We can identify the expression (22) for 7 with the Lagrange's Inversion formula applied at the point z — a with 
G{z) = — \og{z). We obtain 



7 = - log(u') + log(a) . 



(A6) 



We can also compare Lagrange's Inversion formula with the expression (23) for E, taking now G(z) = ^^-j-. We 
obtain 



E 



w+1 a+1 



(A7) 



Eliminating w between the last two equations leads to the desired expression ( 26 ) . We note that the precise form 
of 4>{z) has played no role. This suggest that the formula (26) for the cumulant generating function should be rather 
universal. 
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